# The Impact of Experiences with COVID-19 on the 2020 Presidential Election and Support for Health Reform
# N.M. Kavanagh, A. Menon
# December 6, 2024

# Please direct questions about this script file to nolankavanagh@fas.harvard.edu.

# Clear R environment
rm(list = ls())

# Load packages
library(here)         # Working directory
library(dplyr)        # Analysis tools
library(tidyr)        # Analysis tools
library(stringr)      # Analysis tools
library(lubridate)    # Analysis tools
library(survey)       # Analysis tools
library(arsenal)      # Table tools
library(modelsummary) # Table tools
library(kableExtra)   # Table tools
library(estimatr)     # Modeling tools
library(ggplot2)      # Graphing tools
library(ggforce)      # Graphing tools
library(scales)       # Graphing tools
library(usmap)        # Graphing tools
library(cowplot)      # Graphing tools
library(viridis)      # Graphing tools
library(readstata13)  # Import dataset

##############################################################################
# CES 2020 preparation
##############################################################################

# Read 2020 dataset into R
cces20 <- read.csv("Data/CES 2020/CES20_Common_OUTPUT_vv.csv")
cces20$year <- "2020"

# Reformat interview date
cces20$interview_date <- gsub(" .*", "", cces20$starttime)
cces20$interview_date <- as.Date(cces20$interview_date)

# Age
cces20 <- cces20 %>% mutate(
  age_recode = case_when(
    birthyr %in% c(1990:2002) ~ "18 to 30",
    birthyr %in% c(1980:1989) ~ "31 to 40",
    birthyr %in% c(1970:1979) ~ "41 to 50",
    birthyr %in% c(1960:1969) ~ "51 to 60",
    birthyr %in% c(1900:1959) ~ "61 and older"
  ))

# Gender
cces20 <- cces20 %>% mutate(
  gender_recode = case_when(
    gender == 1 ~ "Male",
    gender == 2 ~ "Female"
  ))

# Race/ethnicity
cces20 <- cces20 %>% mutate(
  race_recode = case_when(
    race == 1        ~ "White",
    race == 2        ~ "Black",
    race == 3        ~ "Hispanic",
    race == 4        ~ "Asian",
    race == 5        ~ "Native American",
    race %in% c(6:8) ~ "Mixed/other"
  ))
cces20$race_recode <- factor(cces20$race_recode, levels = c("White", "Black", "Hispanic", "Asian", "Native American", "Mixed/other"))

# Education
cces20 <- cces20 %>% mutate(
  educ_recode = case_when(
    educ == 1        ~ "Less than high school",
    educ == 2        ~ "High school graduate",
    educ %in% c(3:4) ~ "Some college/2-year degree",
    educ == 5        ~ "College graduate",
    educ == 6        ~ "Postgraduate degree"
  ))
cces20$educ_recode <- factor(cces20$educ_recode, levels = c("Less than high school", "High school graduate", "Some college/2-year degree", "College graduate", "Postgraduate degree"))

# Employment status
cces20 <- cces20 %>% mutate(
  employment = case_when(
    employ %in% c(1:2) ~ "Employed",   # Full and part-time
    employ %in% c(3:4) ~ "Unemployed", # Temporarily and not
    employ %in% c(5)   ~ "Retired",    # As advertised
    employ %in% c(6)   ~ "Disabled",   # Permanently disabled
    employ %in% c(7:9) ~ "Other"       # Stay-at-home, student, other
  ))
cces20$employment <- factor(cces20$employment, levels = c("Employed", "Unemployed", "Retired", "Disabled", "Other"))

# Family income
cces20 <- cces20 %>% mutate(
  income = case_when(
    faminc_new %in% c(1:2)   ~ "$19,999 or less",
    faminc_new %in% c(3:4)   ~ "$20,000 to $39,999",
    faminc_new %in% c(5:6)   ~ "$40,000 to $59,999",
    faminc_new %in% c(7:8)   ~ "$60,000 to $79,999",
    faminc_new %in% c(9)     ~ "$80,000 to $99,999",
    faminc_new %in% c(10:16) ~ "$100,000 or more"
  ))
cces20$income <- factor(cces20$income, levels = c("$19,999 or less", "$20,000 to $39,999", "$40,000 to $59,999", "$60,000 to $79,999", "$80,000 to $99,999", "$100,000 or more"))

# Insurance status
cces20 <- cces20 %>% mutate(
  insurance = case_when (
    healthins_1 == 1 | healthins_3 == 1 ~ "Employer-based",   # Private takes precedence
    healthins_4 == 1                    ~ "Purchased/exchanges",
    healthins_2 == 1                    ~ "Government-based",
    healthins_5 == 1 | healthins_6 == 1 ~ "Unsure or none"
  ))
cces20$insurance <- factor(cces20$insurance, levels = c("Employer-based", "Purchased/exchanges", "Government-based", "Unsure or none"))

# Party identification
cces20 <- cces20 %>% mutate(
  party_id = case_when(
    pid3 == 1        ~ "Democrat",
    pid3 == 2        ~ "Republican",
    pid3 == 3        ~ "Independent",
    pid3 %in% c(4:5) ~ "Other/unsure"
  ))
cces20$party_id <- factor(cces20$party_id, levels = c("Independent", "Democrat", "Republican", "Other/unsure"))

# Strength of partisanship
cces20 <- cces20 %>% mutate(
  party_lean = case_when(
    pid7 %in% c(1,7)     ~ "Strong partisan",
    pid7 %in% c(2:6,8:9) ~ "Weak/non-partisan"
  ))
cces20$party_lean <- factor(cces20$party_lean, levels = c("Weak/non-partisan", "Strong partisan"))

# Know someone with COVID-19
# Main, fully combined measure
cces20 <- cces20 %>% mutate(
  know_covid = case_when(
    CC20_309a_1 == 1 ~ 1, # Self
    CC20_309a_2 == 1 ~ 1, # Family member
    CC20_309a_3 == 1 ~ 1, # Friend
    CC20_309a_4 == 1 ~ 1, # Co-worker
    CC20_309a_5 == 1 ~ 0, # No one
    TRUE             ~ 1  # Anyone else
  ))

# Know someone with COVID-19
# Define as separate persons
cces20 <- cces20 %>% mutate(
  covid_self = case_when(
    CC20_309a_1 == 1 ~ 1, # Self
    TRUE             ~ 0  # Not self
  ))
cces20 <- cces20 %>% mutate(
  covid_family = case_when(
    CC20_309a_2 == 1 ~ 1, # Family member
    TRUE             ~ 0  # Not family
  ))
cces20 <- cces20 %>% mutate(
  covid_friend = case_when(
    CC20_309a_3 == 1 ~ 1, # Friend
    TRUE             ~ 0  # Not friend
  ))
cces20 <- cces20 %>% mutate(
  covid_coworker = case_when(
    CC20_309a_4 == 1 ~ 1, # Co-worker
    TRUE             ~ 0  # Not coworker
  ))
cces20 <- cces20 %>% mutate(
  covid_any_other = case_when(
    CC20_309a_2 == 1 ~ 1, # Family member
    CC20_309a_3 == 1 ~ 1, # Friend
    CC20_309a_4 == 1 ~ 1, # Co-worker
    TRUE             ~ 0  # Not family, friend, coworker
  ))

# Know someone who died of COVID-19
cces20 <- cces20 %>% mutate(
  death_family = case_when(
    CC20_309b_1 == 1 ~ 1, # Family member
    TRUE             ~ 0  # Not family
  ))
cces20 <- cces20 %>% mutate(
  death_friend = case_when(
    CC20_309b_2 == 1 ~ 1, # Friend
    TRUE             ~ 0  # Not friend
  ))
cces20 <- cces20 %>% mutate(
  death_coworker = case_when(
    CC20_309b_3 == 1 ~ 1, # Co-worker
    TRUE             ~ 0  # Not coworker
  ))
cces20 <- cces20 %>% mutate(
  death_any_other = case_when(
    CC20_309b_1 == 1 ~ 1, # Family member
    CC20_309b_2 == 1 ~ 1, # Friend
    CC20_309b_3 == 1 ~ 1, # Co-worker
    TRUE             ~ 0  # Not family, friend, coworker
  ))

# Medicare for All support
cces20 <- cces20 %>% mutate(
  medi_exp = case_when(
    CC20_327a == 1 ~ 1, # Support
    CC20_327a == 2 ~ 0  # Oppose
  ))

# Support for Trump
cces20 <- cces20 %>% mutate(
  trump_2020 = case_when(
    
    # Candidate you prefer?
    CC20_364b == 1        ~ 1,
    CC20_364b %in% c(2:5) ~ 0,
    
    # If already voted, whom?
    CC20_364a == 1        ~ 1,
    CC20_364a %in% c(2:5) ~ 0
  ))

##############################################################################
# CES 2019 data preparation
##############################################################################

# Read 2019 dataset into R
cces19 <- read.dta13("Data/CES 2019/CCES19_Common_OUTPUT.dta")
cces19$year <- "2019"

# Reformat interview date
cces19$interview_date <- gsub(" .*", "", cces19$starttime)
cces19$interview_date <- as.Date(cces19$interview_date)

# Age
cces19 <- cces19 %>% mutate(
  age_recode = case_when(
    birthyr %in% c(1989:2001) ~ "18 to 30",
    birthyr %in% c(1979:1988) ~ "31 to 40",
    birthyr %in% c(1969:1978) ~ "41 to 50",
    birthyr %in% c(1959:1968) ~ "51 to 60",
    birthyr %in% c(1900:1958) ~ "61 and older"
  ))

# Gender
cces19 <- cces19 %>% mutate(
  gender_recode = case_when(
    gender == "Male"   ~ "Male",
    gender == "Female" ~ "Female"
  ))

# Race/ethnicity
cces19 <- cces19 %>% mutate(
  race_recode = case_when(
    race == "White"           ~ "White",
    race == "Black"           ~ "Black",
    race == "Hispanic"        ~ "Hispanic",
    race == "Asian"           ~ "Asian",
    race == "Native American" ~ "Native American",
    race == "Middle Eastern"  ~ "Mixed/other",
    race == "Mixed"           ~ "Mixed/other",
    race == "Other"           ~ "Mixed/other"
  ))
cces19$race_recode <- factor(cces19$race_recode, levels = c("White", "Black", "Hispanic", "Asian", "Native American", "Mixed/other"))

# Education
cces19 <- cces19 %>% mutate(
  educ_recode = case_when(
    educ == "No HS"                ~ "Less than high school",
    educ == "High school graduate" ~ "High school graduate",
    educ == "Some college"         ~ "Some college/2-year degree",
    educ == "2-year"               ~ "Some college/2-year degree",
    educ == "4-year"               ~ "College graduate",
    educ == "Post-grad"            ~ "Postgraduate degree"
  ))
cces19$educ_recode <- factor(cces19$educ_recode, levels = c("Less than high school", "High school graduate", "Some college/2-year degree", "College graduate", "Postgraduate degree"))

# Employment status
cces19 <- cces19 %>% mutate(
  employment = case_when(
    employ %in% c("Full-time", "Part-time")             ~ "Employed",
    employ %in% c("Temporarily laid off", "Unemployed") ~ "Unemployed",
    employ %in% c("Retired")                            ~ "Retired",
    employ %in% c("Permanently disabled")               ~ "Disabled",
    employ %in% c("Homemaker", "Student", "Other")      ~ "Other"
  ))
cces19$employment <- factor(cces19$employment, levels = c("Employed", "Unemployed", "Retired", "Disabled", "Other"))

# Family income
cces19 <- cces19 %>% mutate(
  income = case_when(
    faminc_new %in% c("Less than $10,000", "$10,000 - $19,999")   ~ "$19,999 or less",
    faminc_new %in% c("$20,000 - $29,999", "$30,000 - $39,999")   ~ "$20,000 to $39,999",
    faminc_new %in% c("$40,000 - $49,999", "$50,000 - $59,999")   ~ "$40,000 to $59,999",
    faminc_new %in% c("$60,000 - $69,999", "$70,000 - $79,999")   ~ "$60,000 to $79,999",
    faminc_new %in% c("$80,000 - $99,999")                        ~ "$80,000 to $99,999",
    faminc_new %in% c("$100,000 - $119,999", "$100,000 - $119,999",
                      "$120,000 - $149,999", "$150,000 - $199,999",
                      "$200,000 - $249,999", "$250,000 - $349,999",
                      "$350,000 - $499,999", "$500,000 or more")  ~ "$100,000 or more"
  ))
cces19$income <- factor(cces19$income, levels = c("$19,999 or less", "$20,000 to $39,999", "$40,000 to $59,999", "$60,000 to $79,999", "$80,000 to $99,999", "$100,000 or more"))

# Insurance status
cces19 <- cces19 %>% mutate(
  insurance = case_when (
    healthins_1 == "selected" | healthins_3 == "selected" ~ "Employer-based",
    healthins_4 == "selected"                             ~ "Purchased/exchanges",
    healthins_2 == "selected"                             ~ "Government-based",
    healthins_5 == "selected" | healthins_6 == "selected" ~ "Unsure or none"
  ))
cces19$insurance <- factor(cces19$insurance, levels = c("Employer-based", "Purchased/exchanges", "Government-based", "Unsure or none"))

# Medicare for All support
cces19 <- cces19 %>% mutate(
  medi_exp = case_when(
    CC19_327a == "Support" ~ 1,
    CC19_327a == "Oppose"  ~ 0
  ))

# Party identification
cces19 <- cces19 %>% mutate(
  party_id = case_when(
    pid3 == "Democrat"               ~ "Democrat",
    pid3 == "Republican"             ~ "Republican",
    pid3 == "Independent"            ~ "Independent",
    pid3 %in% c("Other", "Not sure") ~ "Other/unsure"
  ))
cces19$party_id <- factor(cces19$party_id, levels = c("Independent", "Democrat", "Republican", "Other/unsure"))

# Recalled support for Trump in 2016
cces19 <- cces19 %>% mutate(
  trump_2016 = case_when(
    presvote16post == "Donald Trump" ~ 1,
    !is.na(presvote16post)           ~ 0
  ))

# Separate state variable for 2019
cces19$state_alt <- cces19$inputstate

# Reformat FIPS variable
cces19$countyfips <- as.numeric(as.character(cces19$countyfips))

# Select relevant variables
cces19 <- cces19 %>% select(year, interview_date, gender_recode, age_recode, race_recode, educ_recode, employment, income, insurance, medi_exp, party_id, trump_2016, state_alt, countyfips, commonweight)

# Merge datasets
cces_merge <- cces20 %>% full_join(cces19)

##############################################################################
# COVID-19 data preparation
##############################################################################

# Read COVID-19 datasets into R
covid     <- read.csv("Data/COVID-19/us-counties.csv")
covid_nyc <- read.csv("Data/COVID-19/cases-by-day_nyc_counties.csv")

# Reformat NYC dates
covid_nyc$date <- as.Date(covid_nyc$date_of_interest, "%m/%d/%Y")

# Sum daily cases in NYC boroughs
# Returns cumulative cases on each day
covid_nyc[,"bronx"]     <- cumsum(covid_nyc$BX_CASE_COUNT)
covid_nyc[,"manhattan"] <- cumsum(covid_nyc$MN_CASE_COUNT)
covid_nyc[,"kings"]     <- cumsum(covid_nyc$BK_CASE_COUNT)
covid_nyc[,"queens"]    <- cumsum(covid_nyc$QN_CASE_COUNT)
covid_nyc[,"staten"]    <- cumsum(covid_nyc$SI_CASE_COUNT)

# Select matching date ranges
covid     <- subset(covid,     as.Date(covid$date) > as.Date("2020-02-28"))
covid_nyc <- subset(covid_nyc, covid_nyc$date      < as.Date("2021-07-17"))

# Select relevant columns
covid     <- covid     %>% select(date, county, state, fips, cases)
covid_nyc <- covid_nyc %>% select(date, bronx, manhattan, kings, queens, staten)  

# Spread master dataset
covid <- covid %>% spread(date, cases)

# Transpose NYC dataset
covid_nyc <- setNames(data.frame(t(covid_nyc[,-1])), covid_nyc[,1])

# Assign NYC borough information
covid_nyc$county <- c("bronx", "manhattan", "kings", "queens", "staten")
covid_nyc$state  <- "New York"    
covid_nyc$fips   <- c(36005, 36061, 36047, 36081, 36085)

# Merge COVID-19 datasets
covid <- rbind.data.frame(covid, covid_nyc)

# Set NAs to 0, i.e. before first reported cases
covid[is.na(covid)] <- 0

# Read Census dataset (DP05) into R
# Includes county populations for 2015-2019
census <- read.csv("Data/ACSDP5Y2019.DP05_data_with_overlays_2021-07-23T150146.csv")

# Remove extraneous row
census <- census[-grep("id", census$GEO_ID),]

# Reformat FIPS codes
census$GEO_ID <- as.numeric(str_remove(census$GEO_ID, "0500000US"))

# Merge COVID-19 and Census data
covid <- merge(covid, census, by.x = "fips", by.y = "GEO_ID")

# Merge with CES data
cces <- cces_merge %>% merge(covid, by.x = "countyfips", by.y = "fips")

##############################################################################
# Assign COVID-19 cases on interview date
##############################################################################

# Create case variables
cces$cases_on_interview  <- NULL
cces$cases_6_mos_earlier <- NULL
cces$cases_6_mos_later   <- NULL

# Select cases on given dates
for (i in 1:nrow(cces)) {
  
  # Get interview date for given row
  date <- as.Date(cces[i, "interview_date"])
  
  # For respondents in 2019, add one year
  if (cces[i, "year"] == "2019") {
    date <- ymd(as.Date(date) %m+% years(1))
  }
  
  # Find columns of relevant dates
  col1 <- grep(date,                              colnames(cces)) # Date of interview
  col2 <- grep(ymd(as.Date(date) %m-% months(6)), colnames(cces)) # 6 months earlier
  col3 <- grep(ymd(as.Date(date) %m+% months(6)), colnames(cces)) # 6 months earlier
  
  # Save cases on those dates
  cces[i, "cases_on_interview"]  <- cces[i, col1]
  cces[i, "cases_6_mos_earlier"] <- cces[i, col2]
  cces[i, "cases_6_mos_later"]   <- cces[i, col3]
}

# Calculate cases per capita
# Divide by 5 to scale regression coefficients as 5 per 100
cces$DP05_0001E <- as.numeric(as.character(cces$DP05_0001E))
cces$covid_pc_int <- cces$cases_on_interview /cces$DP05_0001E*100/5 # 5 per 100
cces$covid_pc_6be <- cces$cases_6_mos_earlier/cces$DP05_0001E*100/5 # 5 per 100
cces$covid_pc_6af <- cces$cases_6_mos_later  /cces$DP05_0001E*100/5 # 5 per 100

##############################################################################
# Full dataset preparation
##############################################################################

# Complete cases for main models' variables
cces_model <- cces %>%
  filter_at(vars(medi_exp, covid_pc_int, gender_recode, age_recode, race_recode,
                 educ_recode, employment, income, insurance, party_id, commonweight),
            all_vars(!is.na(.)))

# Treat fixed effects as factors
cces_model$inputstate     <- as.factor(cces_model$inputstate)
cces_model$state_alt      <- as.factor(cces_model$state_alt)
cces_model$interview_date <- as.factor(cces_model$interview_date)
cces_model$countyfips     <- as.factor(cces_model$countyfips)

# Lists of covariates for main models
cov_19 <- "+ age_recode + gender_recode + race_recode + educ_recode + employment + income + insurance + party_id + state_alt + interview_date"
cov_20 <- "+ age_recode + gender_recode + race_recode + educ_recode + employment + income + insurance + party_id + inputstate + interview_date"

# Variable names for regression tables
var_names <- 
  c("know_covid"    = "Respondent got COVID-19 or knows someone who did",
    "covid_pc_int"  = "Cumulative county-level COVID-19 cases (5 per 100)",
    "covid_pc_6be"  = "Cumulative county-level COVID-19 cases (5 per 100)",
    "covid_pc_6af"  = "Cumulative county-level COVID-19 cases (5 per 100)",
    "age_recode31 to 40"     = "Age: 31 to 40",
    "age_recode41 to 50"     = "Age: 41 to 50",
    "age_recode51 to 60"     = "Age: 51 to 60",
    "age_recode61 and older" = "Age: 61 and older",
    "gender_recodeMale"      = "Gender: Male",
    "race_recodeBlack"           = "Race/ethnicity: Black",
    "race_recodeHispanic"        = "Race/ethnicity: Hispanic",
    "race_recodeAsian"           = "Race/ethnicity: Asian",
    "race_recodeNative American" = "Race/ethnicity: Native American",
    "race_recodeMixed/other"     = "Race/ethnicity: Mixed/other",
    "educ_recodeHigh school graduate"       = "Education: High school graduate",
    "educ_recodeSome college/2-year degree" = "Education: Some college/2-year degree",
    "educ_recodeCollege graduate"           = "Education: College graduate",
    "educ_recodePostgraduate degree"        = "Education: Postgraduate degree",
    "employmentUnemployed" = "Employment: Unemployed",
    "employmentRetired"    = "Employment: Retired",
    "employmentDisabled"   = "Employment: Disabled",
    "employmentOther"      = "Employment: Other",
    "income$20,000 to $39,999" = "Family income: $20,000 to $39,999",
    "income$40,000 to $59,999" = "Family income: $40,000 to $59,999",
    "income$60,000 to $79,999" = "Family income: $60,000 to $79,999",
    "income$80,000 to $99,999" = "Family income: $80,000 to $99,999",
    "income$100,000 or more"   = "Family income: $100,000 or more",
    "insurancePurchased/exchanges" = "Insurance status: Purchased/exchanges",
    "insuranceGovernment-based"    = "Insurance status: Government-based",
    "insuranceUnsure or none"      = "Insurance status: Unsure or none",
    "party_idDemocrat"     = "Party identification: Democrat",
    "party_idRepublican"   = "Party identification: Republican",
    "party_idIndependent"  = "Party identification: Independent",
    "party_idOther/unsure" = "Party identification: Other/unsure")

# Set working directory for exports
setwd("./Exhibits")

##############################################################################
# Demographic characteristics and descriptives
##############################################################################

# Respondent characteristics: unweighted
summary(tableby(year ~ age_recode + gender_recode + race_recode + educ_recode +
                  employment + income + insurance + party_id,
                cces_model, digits.pct=0), text=T)

# Respondent characteristics: weighted
summary(tableby(year ~ age_recode + gender_recode + race_recode + educ_recode +
                  employment + income + insurance + party_id,
                cces_model, digits.pct=0, weights=commonweight), text=T)

# Weighted individual-level measures in 2020
summary(tableby(~medi_exp + trump_2020 + know_covid + covid_pc_int +
                  covid_self + covid_family + covid_friend + covid_coworker +
                  covid_any_other + death_family + death_friend + death_coworker,
                subset(cces_model, year=="2020"),
                digits=2, weights=commonweight), text=T)

# Weighted Medicare for All support by party in 2020
summary(tableby(party_id ~ medi_exp,
                subset(cces_model, year=="2020"),
                digits=2, weights=commonweight), text=T)

# Weighted Medicare for All support by party in 2019
summary(tableby(party_id ~ medi_exp,
                subset(cces_model, year=="2019"),
                digits=2, weights=commonweight), text=T)

# Range of interview dates
table(cces_model$interview_date)

##############################################################################
# Medicare for All: Main models
##############################################################################

# Main model 1: Known COVID-19 exposure
m4a_main_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Main model 2: Cases per capita
m4a_main_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Main model 3: Both measures
m4a_main_3 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid + covid_pc_int", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1" = m4a_main_1,
                  "Model 2" = m4a_main_2,
                  "Model 3" = m4a_main_3),
             coef_omit   = "state|date|Intercept",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_map    = var_names,
             add_rows    =
               as.data.frame(
                 rbind(cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Full results for specifications reported under Table 1.",
             stars       = FALSE,
             output      = "M4A table for main models.tex") %>%
  row_spec(row=c(96,98), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 3))

##############################################################################
# M4A robustness check: Placebo tests
##############################################################################

# Robustness model: M4A surveyed in 2019, cases after 12 mos.
m4a_plc_1 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_19)),
                       data     = subset(cces_model, year=="2019"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Robustness model: M4A surveyed in 2019, cases after 6 mos.
m4a_plc_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_6be", cov_19)),
                       data     = subset(cces_model, year=="2019"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Robustness model: M4A surveyed in 2019, cases after 18 mos.
m4a_plc_3 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_6af", cov_19)),
                       data     = subset(cces_model, year=="2019"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Extract coefficients and CIs
m4a_CIs <- as.data.frame(rbind(
  cbind("2020\nmain",          m4a_main_2$coefficients[[2]],
        m4a_main_2$conf.low[[2]], m4a_main_2$conf.high[[2]]),
  cbind("2019\n12 mos. after", m4a_plc_1$coefficients[[2]],
        m4a_plc_1$conf.low[[2]], m4a_plc_1$conf.high[[2]]),
  cbind("2019\n6 mos. after",  m4a_plc_2$coefficients[[2]],
        m4a_plc_2$conf.low[[2]], m4a_plc_2$conf.high[[2]]),
  cbind("2019\n18 mos. after", m4a_plc_3$coefficients[[2]],
        m4a_plc_3$conf.low[[2]], m4a_plc_3$conf.high[[2]])
))
colnames(m4a_CIs) <- c("model", "coefficient", "conf_low", "conf_high")
m4a_CIs$model <- factor(m4a_CIs$model, levels = c("2019\n18 mos. after", "2019\n6 mos. after", "2019\n12 mos. after", "2020\nmain"))

# Treat as numeric
m4a_CIs$coefficient <- as.numeric(as.character(m4a_CIs$coefficient))
m4a_CIs$conf_low    <- as.numeric(as.character(m4a_CIs$conf_low))
m4a_CIs$conf_high   <- as.numeric(as.character(m4a_CIs$conf_high))

# Coefficient plot
m4a_coef_plot <- ggplot(m4a_CIs, aes(x=model, y=coefficient,
                                     ymin=conf_low, ymax=conf_high)) +
  geom_point() +
  geom_errorbar(width=0.2) +
  theme_test() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        text = element_text(size = 10, face = "bold"),
        axis.ticks   = element_blank(),
        axis.text.y  = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25)) +
  coord_flip() +
  ylab("\nPercentage-point difference in\nM4A support given 5 per 100 increase\nin COVID-19 cases per capita") +
  geom_hline(yintercept=0, color="black") +
  scale_y_continuous(labels=scales::percent_format(accuracy=1),
                     breaks=c(-0.15,-0.1,-0.05,0,0.05,0.1,0.15),
                     limits=c(-0.15,0.15))

# Print figure
ggsave(plot=m4a_coef_plot, file="M4A coefficient plot, weighted.png",
       width=2.7, height=3.75, units='in', dpi=600)

# Compile results into table
modelsummary(list("Line 1" = m4a_main_2,
                  "Line 2" = m4a_plc_1,
                  "Line 3" = m4a_plc_2,
                  "Line 4" = m4a_plc_3),
             coef_omit   = "state|date|Intercept",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_map    = var_names,
             title       = "Full results for specifications reported under Figure 3A.",
             add_rows    =
               as.data.frame(
                 rbind(cbind("State fixed effects",          "Yes", "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes", "Yes"))),
             stars       = FALSE,
             output      = "M4A table for placebo tests.tex") %>%
  row_spec(row=c(93,95), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "2020 CES" = 1, "2019 CES" = 3)) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 4))

##############################################################################
# M4A robustness check: Separate individual exposures
##############################################################################

# Main model 1a: Respondent had COVID-19
m4a_alt_1a <- lm_robust(as.formula(paste0("medi_exp ~ covid_self", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Main model 1b: Any contact had COVID-19
m4a_alt_1b <- lm_robust(as.formula(paste0("medi_exp ~ covid_any_other", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Main model 1c: Self vs. any contact had COVID-19
m4a_alt_1c <- lm_robust(as.formula(paste0("medi_exp ~ covid_self + covid_any_other", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Main model 1d: Separate terms for all persons
m4a_alt_1d <- lm_robust(as.formula(paste0("medi_exp ~ covid_self + covid_family + covid_friend + covid_coworker", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1a" = m4a_alt_1a,
                  "Model 1b" = m4a_alt_1b,
                  "Model 1c" = m4a_alt_1c,
                  "Model 1d" = m4a_alt_1d),
             coef_omit   = "^(?!covid)",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = c("covid_self"      = "Respondent got COVID-19 themselves",
                             "covid_any_other" = "Family, friend, or coworker got COVID-19",
                             "covid_family"    = "Family member got COVID-19",
                             "covid_friend"    = "Friend got COVID-19",
                             "covid_coworker"  = "Coworker got COVID-19"),
             add_rows    =
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes", "Yes"))),
             title       = "Medicare for All Model 1 with separate individual-level exposures.",
             stars       = FALSE,
             output      = "M4A table for separate exposures.tex") %>%
  row_spec(row=c(18,21), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 4))

##############################################################################
# M4A robustness check: Individual-level deaths
##############################################################################

# Main model 1w: Death of family member
m4a_alt_1w <- lm_robust(as.formula(paste0("medi_exp ~ death_family", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Main model 1x: Death of friend
m4a_alt_1x <- lm_robust(as.formula(paste0("medi_exp ~ death_friend", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Main model 1y: Death of coworker
m4a_alt_1y <- lm_robust(as.formula(paste0("medi_exp ~ death_coworker", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Main model 1y: Death of any person (separately)
m4a_alt_1z <- lm_robust(as.formula(paste0("medi_exp ~ death_family + death_friend + death_coworker", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1w" = m4a_alt_1w,
                  "Model 1x" = m4a_alt_1x,
                  "Model 1y" = m4a_alt_1y,
                  "Model 1z" = m4a_alt_1z),
             coef_omit   = "^(?!death)",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = c("death_family"    = "Death of family member due to COVID-19",
                             "death_friend"    = "Death of friend due to COVID-19",
                             "death_coworker"  = "Death of coworker due to COVID-19"),
             add_rows    =
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes", "Yes"))),
             title       = "Medicare for All Model 1 with deaths of individual-level contacts.",
             stars       = FALSE,
             output      = "M4A table for deaths.tex") %>%
  row_spec(row=c(12,15), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 4))

##############################################################################
# M4A robustness check: Without weights
##############################################################################

# Model 1: Known COVID-19 exposure
m4a_mn_uw_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_20)),
                         data     = subset(cces_model, year=="2020"),
                         clusters = countyfips,
                         se_type  = "stata")

# Model 2: Cases per capita
m4a_mn_uw_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_20)),
                         data     = subset(cces_model, year=="2020"),
                         clusters = countyfips,
                         se_type  = "stata")

# Model 3: Both measures
m4a_mn_uw_3 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid + covid_pc_int", cov_20)),
                         data     = subset(cces_model, year=="2020"),
                         clusters = countyfips,
                         se_type  = "stata")

# Robustness model: M4A surveyed in 2019, cases after 12 mos.
m4a_pl_uw_1 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_19)),
                         data     = subset(cces_model, year=="2019"),
                         clusters = countyfips,
                         se_type  = "stata")

# Robustness model: M4A surveyed in 2019, cases after 6 mos.
m4a_pl_uw_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_6be", cov_19)),
                         data     = subset(cces_model, year=="2019"),
                         clusters = countyfips,
                         se_type  = "stata")

# Robustness model: M4A surveyed in 2019, cases after 18 mos.
m4a_pl_uw_3 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_6af", cov_19)),
                         data     = subset(cces_model, year=="2019"),
                         clusters = countyfips,
                         se_type  = "stata")

# Extract CIs
m4a_uw_CIs <- as.data.frame(rbind(
  cbind("2020\nmain",          m4a_mn_uw_2$coefficients[[2]],
        m4a_mn_uw_2$conf.low[[2]],  m4a_mn_uw_2$conf.high[[2]]),
  cbind("2019\n12 mos. after", m4a_pl_uw_1$coefficients[[2]],
        m4a_pl_uw_1$conf.low[[2]], m4a_pl_uw_1$conf.high[[2]]),
  cbind("2019\n6 mos. after",  m4a_pl_uw_2$coefficients[[2]],
        m4a_pl_uw_2$conf.low[[2]], m4a_pl_uw_2$conf.high[[2]]),
  cbind("2019\n18 mos. after", m4a_pl_uw_3$coefficients[[2]],
        m4a_pl_uw_3$conf.low[[2]], m4a_pl_uw_3$conf.high[[2]])
))
colnames(m4a_uw_CIs) <- c("model", "coefficient", "conf_low", "conf_high")
m4a_uw_CIs$model <- factor(m4a_uw_CIs$model, levels = c("2019\n18 mos. after", "2019\n6 mos. after", "2019\n12 mos. after", "2020\nmain"))

# Treat as numeric
m4a_uw_CIs$coefficient <- as.numeric(as.character(m4a_uw_CIs$coefficient))
m4a_uw_CIs$conf_low    <- as.numeric(as.character(m4a_uw_CIs$conf_low))
m4a_uw_CIs$conf_high   <- as.numeric(as.character(m4a_uw_CIs$conf_high))

# Coefficient plot
m4a_uw_coef_plot <- ggplot(m4a_uw_CIs, aes(x=model, y=coefficient,
                                           ymin=conf_low, ymax=conf_high)) +
  geom_point() +
  geom_errorbar(width=0.2) +
  theme_test() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        text = element_text(size = 10, face = "bold"),
        axis.ticks   = element_blank(),
        axis.text.y  = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25)) +
  coord_flip() +
  ylab("\nPercentage-point difference in\nM4A support given 5 per 100 increase\nin COVID-19 cases per capita") +
  geom_hline(yintercept=0, color="black") +
  scale_y_continuous(labels=scales::percent_format(accuracy=1),
                     breaks=c(-0.15,-0.1,-0.05,0,0.05,0.1,0.15),
                     limits=c(-0.15,0.15))

# Print figure
ggsave(plot=m4a_uw_coef_plot, file="M4A coefficient plot, unweighted.png",
       width=2.7, height=3.75, units='in', dpi=600)

# Compile results into table
modelsummary(list("Model 1" = m4a_mn_uw_1,
                  "Model 2" = m4a_mn_uw_2,
                  "Model 3" = m4a_mn_uw_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = 
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "No",  "No",  "No"))),
             title       = "Medicare for All models without survey weights.",
             stars       = FALSE,
             output      = "M4A table for unweighted.tex") %>%
  row_spec(row=c(9,12), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 3))

##############################################################################
# M4A robustness check: Logistic regression
##############################################################################

# Specify design with weights
design_2020 <- svydesign(ids=~countyfips, weights=~commonweight,
                         data=subset(cces_model, year=="2020"))

# Model 1: Known COVID-19 exposure
m4a_glm_1 <- svyglm(as.formula(paste0("medi_exp ~ know_covid", cov_20)),
                    design = design_2020,
                    family = "binomial")

# Model 2: Cases per capita
m4a_glm_2 <- svyglm(as.formula(paste0("medi_exp ~ covid_pc_int", cov_20)),
                    design = design_2020,
                    family = "binomial")

# Model 3: Both measures
m4a_glm_3 <- svyglm(as.formula(paste0("medi_exp ~ know_covid + covid_pc_int", cov_20)),
                    design = design_2020,
                    family = "binomial")

# Compile results into table
modelsummary(list("Model 1" = m4a_glm_1,
                  "Model 2" = m4a_glm_2,
                  "Model 3" = m4a_glm_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std|Adj",
             coef_rename = var_names,
             add_rows    = 
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Medicare for All models using binomial logistic regression.",
             stars       = FALSE,
             output      = "M4A table for logistic.tex") %>%
  row_spec(row=c(8,11), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 3))

##############################################################################
# M4A robustness check: County fixed effects
##############################################################################

# Model 1: Known COVID-19 exposure
m4a_cfe_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid + countyfips", cov_20)),
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 2: Cases per capita
m4a_cfe_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int + countyfips", cov_20)),
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 3: Both measures
m4a_cfe_3 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid + covid_pc_int + countyfips", cov_20)),
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1" = m4a_cfe_1,
                  "Model 2" = m4a_cfe_2,
                  "Model 3" = m4a_cfe_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = 
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("County fixed effects",         "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Medicare for All models with county fixed effects.",
             stars       = FALSE,
             output      = "M4A table for county FEs.tex") %>%
  row_spec(row=c(9,12), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 3))

##############################################################################
# M4A robustness check: Sequential addition of controls
##############################################################################

cov_fes <- "+ inputstate + interview_date"
cov_dem <- "+ age_recode + gender_recode + race_recode + educ_recode"
cov_ses <- "+ employment + income + insurance"

# Model 1: Only demographics
m4a_seq_1_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_fes, cov_dem)),
                         data     = subset(cces_model, year=="2020"),
                         weights  = commonweight,
                         clusters = countyfips,
                         se_type  = "stata")

# Model 1: Demographics + SES
m4a_seq_1_2 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_fes, cov_dem, cov_ses)),
                         data     = subset(cces_model, year=="2020"),
                         weights  = commonweight,
                         clusters = countyfips,
                         se_type  = "stata")

# Model 1: Full model
m4a_seq_1_3 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_20)),
                         data     = subset(cces_model, year=="2020"),
                         weights  = commonweight,
                         clusters = countyfips,
                         se_type  = "stata")

# Additional rows for table
cov_sequential <- as.data.frame(
  rbind(cbind("Demographic controls",         "Yes", "Yes", "Yes"),
        cbind("Socioeconomic controls",       "No",  "Yes", "Yes"),
        cbind("Party identification",         "No",  "No",  "Yes"),
        cbind("State fixed effects",          "Yes", "Yes", "Yes"),
        cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
        cbind("Clustered standard errors", "County", "County", "County"),
        cbind("Survey weights",               "Yes", "Yes", "Yes")))

# Compile results into table
modelsummary(list("Demo."       = m4a_seq_1_1,
                  "Demo. + SES" = m4a_seq_1_2,
                  "Full"        = m4a_seq_1_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = cov_sequential,
             title       = "Medicare for All Model 1 with sequential addition of controls.",
             stars       = FALSE,
             output      = "M4A table for sequential 1.tex") %>%
  row_spec(row=c(6,9,11), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 3))

# Model 2: Only demographics
m4a_seq_2_1 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_fes, cov_dem)),
                         data     = subset(cces_model, year=="2020"),
                         weights  = commonweight,
                         clusters = countyfips,
                         se_type  = "stata")

# Model 2: Demographics + SES
m4a_seq_2_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_fes, cov_dem, cov_ses)),
                         data     = subset(cces_model, year=="2020"),
                         weights  = commonweight,
                         clusters = countyfips,
                         se_type  = "stata")

# Model 2: Full model
m4a_seq_2_3 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_20)),
                         data     = subset(cces_model, year=="2020"),
                         weights  = commonweight,
                         clusters = countyfips,
                         se_type  = "stata")

# Compile results into table
modelsummary(list("Demo."       = m4a_seq_2_1,
                  "Demo. + SES" = m4a_seq_2_2,
                  "Full"        = m4a_seq_2_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = cov_sequential,
             title       = "Medicare for All Model 2 with sequential addition of controls.",
             stars       = FALSE,
             output      = "M4A table for sequential 2.tex") %>%
  row_spec(row=c(6,9,11), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 3))

# Model 3: Only demographics
m4a_seq_3_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid + covid_pc_int", cov_fes, cov_dem)),
                         data     = subset(cces_model, year=="2020"),
                         weights  = commonweight,
                         clusters = countyfips,
                         se_type  = "stata")

# Model 3: Demographics + SES
m4a_seq_3_2 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid + covid_pc_int", cov_fes, cov_dem, cov_ses)),
                         data     = subset(cces_model, year=="2020"),
                         weights  = commonweight,
                         clusters = countyfips,
                         se_type  = "stata")

# Model 3: Full model
m4a_seq_3_3 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid + covid_pc_int", cov_20)),
                         data     = subset(cces_model, year=="2020"),
                         weights  = commonweight,
                         clusters = countyfips,
                         se_type  = "stata")

# Compile results into table
modelsummary(list("Demo."       = m4a_seq_3_1,
                  "Demo. + SES" = m4a_seq_3_2,
                  "Full"        = m4a_seq_3_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = cov_sequential,
             title       = "Medicare for All Model 3 with sequential addition of controls.",
             stars       = FALSE,
             output      = "M4A table for sequential 3.tex") %>%
  row_spec(row=c(9,12,14), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 3))

##############################################################################
# M4A subgroup analyses: Partisan identification
##############################################################################

# Covariates minus party ID
cov_pid <- gsub("+ party_id", "", cov_20)

# Model 1: Democrats only
m4a_dem_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_pid)),
                       data     = subset(cces_model, year=="2020" & party_id=="Democrat"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 2: Democrats only
m4a_dem_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_pid)),
                       data     = subset(cces_model, year=="2020" & party_id=="Democrat"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 1: Republicans only
m4a_rep_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_pid)),
                       data     = subset(cces_model, year=="2020" & party_id=="Republican"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 2: Republicans only
m4a_rep_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_pid)),
                       data     = subset(cces_model, year=="2020" & party_id=="Republican"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 1: Independents only
m4a_ind_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_pid)),
                       data     = subset(cces_model, year=="2020" & party_id=="Independent"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 2: Independents only
m4a_ind_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_pid)),
                       data     = subset(cces_model, year=="2020" & party_id=="Independent"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 1: Parties fully interacted
m4a_pid_1 <- lm_robust(medi_exp ~ know_covid*party_id + age_recode*party_id + gender_recode*party_id + race_recode*party_id + educ_recode*party_id + employment*party_id + income*party_id + insurance*party_id + inputstate*party_id + interview_date*party_id,
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 2: Parties fully interacted
m4a_pid_2 <- lm_robust(medi_exp ~ covid_pc_int*party_id + age_recode*party_id + gender_recode*party_id + race_recode*party_id + educ_recode*party_id + employment*party_id + income*party_id + insurance*party_id + inputstate*party_id + interview_date*party_id,
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 3: Parties fully interacted
m4a_pid_3 <- lm_robust(medi_exp ~ know_covid*party_id + covid_pc_int*party_id + age_recode*party_id + gender_recode*party_id + race_recode*party_id + educ_recode*party_id + employment*party_id + income*party_id + insurance*party_id + inputstate*party_id + interview_date*party_id,
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1" = m4a_pid_1,
                  "Model 2" = m4a_pid_2,
                  "Model 3" = m4a_pid_3),
             coef_omit   = "^(?!.*know_covid|.*covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = c("know_covid" =
                               "Respondent got COVID-19 or knows someone who did",
                             "know_covid:party_idRepublican" =
                               "Got COVID-19 or knows someone × Republican",
                             "party_idRepublican:know_covid" =
                               "Got COVID-19 or knows someone × Republican",
                             "know_covid:party_idDemocrat" =
                               "Got COVID-19 or knows someone × Democrat",
                             "party_idDemocrat:know_covid" =
                               "Got COVID-19 or knows someone × Democrat",
                             "know_covid:party_idOther/unsure" =
                               "Got COVID-19 or knows someone × Other/unsure",
                             "party_idOther/unsure:know_covid" =
                               "Got COVID-19 or knows someone × Other/unsure",
                             "covid_pc_int" =
                               "Cumulative county-level COVID-19 cases (5 per 100)",
                             "covid_pc_int:party_idRepublican" = 
                               "County-level cases (5 per 100) × Republican",
                             "party_idRepublican:covid_pc_int" = 
                               "County-level cases (5 per 100) × Republican",
                             "covid_pc_int:party_idDemocrat" = 
                               "County-level cases (5 per 100) × Democrat",
                             "party_idDemocrat:covid_pc_int" = 
                               "County-level cases (5 per 100) × Democrat",
                             "covid_pc_int:party_idOther/unsure" = 
                               "County-level cases (5 per 100) × Other/unsure",
                             "party_idOther/unsure:covid_pc_int" = 
                               "County-level cases (5 per 100) × Other/unsure"),
             add_rows    =
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Fully interacted Medicare for All models for party identification.",
             stars       = FALSE,
             output      = "M4A table for parties interacted.tex") %>%
  row_spec(row=c(27,30), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 3))

##############################################################################
# M4A subgroup analyses: Strength of partisanship
##############################################################################

# Model 1: Strong partisans only
m4a_str_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_20)),
                       data     = subset(cces_model, year=="2020" & party_lean=="Strong partisan"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 2: Strong partisans only
m4a_str_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_20)),
                       data     = subset(cces_model, year=="2020" & party_lean=="Strong partisan"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 1: Weak/non-partisans only
m4a_weak_1 <- lm_robust(as.formula(paste0("medi_exp ~ know_covid", cov_20)),
                        data     = subset(cces_model, year=="2020" & party_lean=="Weak/non-partisan"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 2: Weak/non-partisans only
m4a_weak_2 <- lm_robust(as.formula(paste0("medi_exp ~ covid_pc_int", cov_20)),
                        data     = subset(cces_model, year=="2020" & party_lean=="Weak/non-partisan"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 1: Partisan strength fully interacted
m4a_lean_1 <- lm_robust(medi_exp ~ know_covid*party_lean + age_recode*party_lean + gender_recode*party_lean + race_recode*party_lean + educ_recode*party_lean + employment*party_lean + income*party_lean + insurance*party_lean + party_id*party_lean + inputstate*party_lean + interview_date*party_lean,
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 2: Partisan strength fully interacted
m4a_lean_2 <- lm_robust(medi_exp ~ covid_pc_int*party_lean + age_recode*party_lean + gender_recode*party_lean + race_recode*party_lean + educ_recode*party_lean + employment*party_lean + income*party_lean + insurance*party_lean + party_id*party_lean + inputstate*party_lean + interview_date*party_lean,
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 3: Partisan strength fully interacted
m4a_lean_3 <- lm_robust(medi_exp ~ know_covid*party_lean + covid_pc_int*party_lean + age_recode*party_lean + gender_recode*party_lean + race_recode*party_lean + educ_recode*party_lean + employment*party_lean + income*party_lean + insurance*party_lean + party_id*party_lean + inputstate*party_lean + interview_date*party_lean,
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1" = m4a_lean_1,
                  "Model 2" = m4a_lean_2,
                  "Model 3" = m4a_lean_3),
             coef_omit   = "^(?!.*know_covid|.*covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = c("know_covid" =
                               "Respondent got COVID-19 or knows someone who did",
                             "know_covid:party_leanStrong partisan" =
                               "Got COVID-19 or knows someone × Strong partisan",
                             "party_leanStrong partisan:know_covid" =
                               "Got COVID-19 or knows someone × Strong partisan",
                             "covid_pc_int" =
                               "Cumulative county-level COVID-19 cases (5 per 100)",
                             "covid_pc_int:party_leanStrong partisan" =
                               "County-level cases (5 per 100) × Strong partisan",
                             "party_leanStrong partisan:covid_pc_int" =
                               "County-level cases (5 per 100) × Strong partisan"),
             add_rows    =
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Fully interacted Medicare for All models for partisan strength.",
             stars       = FALSE,
             output      = "M4A table for strength interacted.tex") %>%
  row_spec(row=c(15,18), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Medicare for All support" = 3))

##############################################################################
# Trump support: Main models
##############################################################################

# Model 1: Known COVID-19 exposure
tr_main_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_20)),
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 2: Cases per capita
tr_main_2 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_20)),
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 3: Both measures
tr_main_3 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid + covid_pc_int", cov_20)),
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1" = tr_main_1,
                  "Model 2" = tr_main_2,
                  "Model 3" = tr_main_3),
             coef_omit   = "^(?!.*know_covid|.*covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    =
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Main models for Trump support.",
             stars       = FALSE,
             output      = "Trump table for main models.tex")  %>%
  row_spec(row=c(9,12), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Support for Trump" = 3))

##############################################################################
# Trump support: Placebo tests
##############################################################################

# Robustness model: Trump support surveyed in 2019, cases after 12 mos.
tr_plc_1 <- lm_robust(as.formula(paste0("trump_2016 ~ covid_pc_int", cov_19)),
                      data     = subset(cces_model, year=="2019"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Robustness model: Trump support surveyed in 2019, cases after 6 mos.
tr_plc_2 <- lm_robust(as.formula(paste0("trump_2016 ~ covid_pc_6be", cov_19)),
                      data     = subset(cces_model, year=="2019"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Robustness model: Trump support surveyed in 2019, cases after 18 mos.
tr_plc_3 <- lm_robust(as.formula(paste0("trump_2016 ~ covid_pc_6af", cov_19)),
                      data     = subset(cces_model, year=="2019"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Extract coefficients and CIs
tr_CIs <- as.data.frame(rbind(
  cbind("2020\nmain",          tr_main_2$coefficients[[2]],
        tr_main_2$conf.low[[2]], tr_main_2$conf.high[[2]]),
  cbind("2019\n12 mos. after", tr_plc_1$coefficients[[2]],
        tr_plc_1$conf.low[[2]], tr_plc_1$conf.high[[2]]),
  cbind("2019\n6 mos. after",  tr_plc_2$coefficients[[2]],
        tr_plc_2$conf.low[[2]], tr_plc_2$conf.high[[2]]),
  cbind("2019\n18 mos. after", tr_plc_3$coefficients[[2]],
        tr_plc_3$conf.low[[2]], tr_plc_3$conf.high[[2]])
))
colnames(tr_CIs) <- c("model", "coefficient", "conf_low", "conf_high")
tr_CIs$model <- factor(tr_CIs$model, levels = c("2019\n18 mos. after", "2019\n6 mos. after", "2019\n12 mos. after", "2020\nmain"))

# Treat as numeric
tr_CIs$coefficient <- as.numeric(as.character(tr_CIs$coefficient))
tr_CIs$conf_low    <- as.numeric(as.character(tr_CIs$conf_low))
tr_CIs$conf_high   <- as.numeric(as.character(tr_CIs$conf_high))

# Coefficient plot
coef_plot_trump <- ggplot(tr_CIs, aes(x=model, y=coefficient,
                                      ymin=conf_low, ymax=conf_high)) +
  geom_point() +
  geom_errorbar(width=0.2) +
  theme_test() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        text = element_text(size = 10, face = "bold"),
        axis.ticks   = element_blank(),
        axis.text.y  = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25)) +
  coord_flip() +
  ylab("\nPercentage-point difference in\nTrump support given 5 per 100 increase\nin COVID-19 cases per capita") +
  geom_hline(yintercept=0, color="black") +
  scale_y_continuous(labels=scales::percent_format(accuracy=1),
                     breaks=c(-0.15,-0.1,-0.05,0,0.05,0.1,0.15),
                     limits=c(-0.15,0.15))

# Print figure
ggsave(plot=coef_plot_trump, file="Trump coefficient plot, weighted.png",
       width=2.7, height=3.75, units='in', dpi=600)

# Compile results into table
modelsummary(list("Line 1" = tr_main_2,
                  "Line 2" = tr_plc_1,
                  "Line 3" = tr_plc_2,
                  "Line 4" = tr_plc_3),
             coef_omit   = "state|date|Intercept",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_map    = var_names,
             title       = "Full results for specifications reported under Figure 3B.",
             add_rows    =
               as.data.frame(
                 rbind(cbind("State fixed effects",          "Yes", "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes", "Yes"))),
             stars       = FALSE,
             output      = "Trump table for placebo tests.tex") %>%
  row_spec(row=c(93,95), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "2020 CES" = 1, "2019 CES" = 3)) %>%
  add_header_above(c(" " = 1, "Trump support" = 4))

##############################################################################
# Trump robustness check: Without weights
##############################################################################

# Model 1: Known COVID-19 exposure
tr_mn_uw_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        clusters = countyfips,
                        se_type  = "stata")

# Model 2: Cases per capita
tr_mn_uw_2 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        clusters = countyfips,
                        se_type  = "stata")

# Model 3: Both measures
tr_mn_uw_3 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid + covid_pc_int", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        clusters = countyfips,
                        se_type  = "stata")

# Robustness model: Trump surveyed in 2019, cases after 12 mos.
tr_pl_uw_1 <- lm_robust(as.formula(paste0("trump_2016 ~ covid_pc_int", cov_19)),
                        data     = subset(cces_model, year=="2019"),
                        clusters = countyfips,
                        se_type  = "stata")

# Robustness model: Trump surveyed in 2019, cases after 6 mos.
tr_pl_uw_2 <- lm_robust(as.formula(paste0("trump_2016 ~ covid_pc_6be", cov_19)),
                        data     = subset(cces_model, year=="2019"),
                        clusters = countyfips,
                        se_type  = "stata")

# Robustness model: Trump surveyed in 2019, cases after 18 mos.
tr_pl_uw_3 <- lm_robust(as.formula(paste0("trump_2016 ~ covid_pc_6af", cov_19)),
                        data     = subset(cces_model, year=="2019"),
                        clusters = countyfips,
                        se_type  = "stata")

# Extract CIs
tr_uw_CIs <- as.data.frame(rbind(
  cbind("2020\nmain",          tr_mn_uw_2$coefficients[[2]],
        tr_mn_uw_2$conf.low[[2]], tr_mn_uw_2$conf.high[[2]]),
  cbind("2019\n12 mos. after", tr_pl_uw_1$coefficients[[2]],
        tr_pl_uw_1$conf.low[[2]], tr_pl_uw_1$conf.high[[2]]),
  cbind("2019\n6 mos. after",  tr_pl_uw_2$coefficients[[2]],
        tr_pl_uw_2$conf.low[[2]], tr_pl_uw_2$conf.high[[2]]),
  cbind("2019\n18 mos. after", tr_pl_uw_3$coefficients[[2]],
        tr_pl_uw_3$conf.low[[2]], tr_pl_uw_3$conf.high[[2]])
))
colnames(tr_uw_CIs) <- c("model", "coefficient", "conf_low", "conf_high")
tr_uw_CIs$model <- factor(tr_uw_CIs$model, levels = c("2019\n18 mos. after", "2019\n6 mos. after", "2019\n12 mos. after", "2020\nmain"))

# Treat as numeric
tr_uw_CIs$coefficient <- as.numeric(as.character(tr_uw_CIs$coefficient))
tr_uw_CIs$conf_low    <- as.numeric(as.character(tr_uw_CIs$conf_low))
tr_uw_CIs$conf_high   <- as.numeric(as.character(tr_uw_CIs$conf_high))

# Coefficient plot
tr_uw_coef_plot <- ggplot(tr_uw_CIs, aes(x=model, y=coefficient,
                                         ymin=conf_low, ymax=conf_high)) +
  geom_point() +
  geom_errorbar(width=0.2) +
  theme_test() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        text = element_text(size = 10, face = "bold"),
        axis.ticks   = element_blank(),
        axis.text.y  = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25)) +
  coord_flip() +
  ylab("\nPercentage-point difference in\nTrump support given 5 per 100 increase\nin COVID-19 cases per capita") +
  geom_hline(yintercept=0, color="black") +
  scale_y_continuous(labels=scales::percent_format(accuracy=1),
                     breaks=c(-0.15,-0.1,-0.05,0,0.05,0.1,0.15),
                     limits=c(-0.15,0.15))

# Print figure
ggsave(plot=tr_uw_coef_plot, file="Trump coefficient plot, unweighted.png",
       width=2.7, height=3.75, units='in', dpi=600)

# Compile results into table
modelsummary(list("Model 1" = tr_mn_uw_1,
                  "Model 2" = tr_mn_uw_2,
                  "Model 3" = tr_mn_uw_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = 
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "No",  "No",  "No"))),
             title       = "Trump models without survey weights.",
             stars       = FALSE,
             output      = "Trump table for unweighted.tex") %>%
  row_spec(row=c(9,12), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Trump support" = 3))

##############################################################################
# Trump robustness check: Logistic regression
##############################################################################

# Model 1: Known COVID-19 exposure
tr_glm_1 <- svyglm(as.formula(paste0("trump_2020 ~ know_covid", cov_20)),
                   design = design_2020,
                   family = "binomial")

# Model 2: Cases per capita
tr_glm_2 <- svyglm(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_20)),
                   design = design_2020,
                   family = "binomial")

# Model 3: Both measures
tr_glm_3 <- svyglm(as.formula(paste0("trump_2020 ~ know_covid + covid_pc_int", cov_20)),
                   design = design_2020,
                   family = "binomial")

# Compile results into table
modelsummary(list("Model 1" = tr_glm_1,
                  "Model 2" = tr_glm_2,
                  "Model 3" = tr_glm_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std|Adj",
             coef_rename = var_names,
             add_rows    = 
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Trump models using binomial logistic regression.",
             stars       = FALSE,
             output      = "Trump table for logistic.tex") %>%
  row_spec(row=c(8,11), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Trump support" = 3))

##############################################################################
# Trump robustness check: County fixed effects
##############################################################################

# Model 1: Known COVID-19 exposure
tr_cfe_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid + countyfips", cov_20)),
                      data     = subset(cces_model, year=="2020"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 2: Cases per capita
tr_cfe_2 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int + countyfips", cov_20)),
                      data     = subset(cces_model, year=="2020"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 3: Both measures
tr_cfe_3 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid + covid_pc_int + countyfips", cov_20)),
                      data     = subset(cces_model, year=="2020"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1" = tr_cfe_1,
                  "Model 2" = tr_cfe_2,
                  "Model 3" = tr_cfe_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = 
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("County fixed effects",         "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Trump models with county fixed effects.",
             stars       = FALSE,
             output      = "Trump table for county FEs.tex") %>%
  row_spec(row=c(9,12), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Trump support" = 3))

##############################################################################
# Trump robustness check: Sequential addition of controls
##############################################################################

cov_fes <- "+ inputstate + interview_date"
cov_dem <- "+ age_recode + gender_recode + race_recode + educ_recode"
cov_ses <- "+ employment + income + insurance"

# Model 1: Only demographics
tr_seq_1_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_fes, cov_dem)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 1: Demographics + SES
tr_seq_1_2 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_fes, cov_dem, cov_ses)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 1: Full model
tr_seq_1_3 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Additional rows for table
cov_sequential <- as.data.frame(
  rbind(cbind("Demographic controls",         "Yes", "Yes", "Yes"),
        cbind("Socioeconomic controls",       "No",  "Yes", "Yes"),
        cbind("Party identification",         "No",  "No",  "Yes"),
        cbind("State fixed effects",          "Yes", "Yes", "Yes"),
        cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
        cbind("Clustered standard errors", "County", "County", "County"),
        cbind("Survey weights",               "Yes", "Yes", "Yes")))

# Compile results into table
modelsummary(list("Demo."       = tr_seq_1_1,
                  "Demo. + SES" = tr_seq_1_2,
                  "Full"        = tr_seq_1_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = cov_sequential,
             title       = "Trump Model 1 with sequential addition of controls.",
             stars       = FALSE,
             output      = "Trump table for sequential 1.tex") %>%
  row_spec(row=c(6,9,11), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Trump support" = 3))

# Model 2: Only demographics
tr_seq_2_1 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_fes, cov_dem)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 2: Demographics + SES
tr_seq_2_2 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_fes, cov_dem, cov_ses)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 2: Full model
tr_seq_2_3 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Compile results into table
modelsummary(list("Demo."       = tr_seq_2_1,
                  "Demo. + SES" = tr_seq_2_2,
                  "Full"        = tr_seq_2_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = cov_sequential,
             title       = "Trump Model 2 with sequential addition of controls.",
             stars       = FALSE,
             output      = "Trump table for sequential 2.tex") %>%
  row_spec(row=c(6,9,11), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Trump support" = 3))

# Model 3: Only demographics
tr_seq_3_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid + covid_pc_int", cov_fes, cov_dem)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 3: Demographics + SES
tr_seq_3_2 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid + covid_pc_int", cov_fes, cov_dem, cov_ses)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Model 3: Full model
tr_seq_3_3 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid + covid_pc_int", cov_20)),
                        data     = subset(cces_model, year=="2020"),
                        weights  = commonweight,
                        clusters = countyfips,
                        se_type  = "stata")

# Compile results into table
modelsummary(list("Demo."       = tr_seq_3_1,
                  "Demo. + SES" = tr_seq_3_2,
                  "Full"        = tr_seq_3_3),
             coef_omit   = "^(?!know_covid|covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("std.error", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = var_names,
             add_rows    = cov_sequential,
             title       = "Trump Model 3 with sequential addition of controls.",
             stars       = FALSE,
             output      = "Trump table for sequential 3.tex") %>%
  row_spec(row=c(9,12,14), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Trump support" = 3))

##############################################################################
# Trump subgroup analyses: Partisan identification
##############################################################################

# Covariates minus party ID
cov_pid <- gsub("+ party_id", "", cov_20)

# Model 1: Democrats only
tr_dem_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_pid)),
                      data     = subset(cces_model, year=="2020" & party_id=="Democrat"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 2: Democrats only
tr_dem_2 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_pid)),
                      data     = subset(cces_model, year=="2020" & party_id=="Democrat"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 1: Republicans only
tr_rep_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_pid)),
                      data     = subset(cces_model, year=="2020" & party_id=="Republican"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 2: Republicans only
tr_rep_2 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_pid)),
                      data     = subset(cces_model, year=="2020" & party_id=="Republican"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 1: Independents only
tr_ind_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_pid)),
                      data     = subset(cces_model, year=="2020" & party_id=="Independent"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 2: Independents only
tr_ind_2 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_pid)),
                      data     = subset(cces_model, year=="2020" & party_id=="Independent"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 1: Parties fully interacted
tr_pid_1 <- lm_robust(trump_2020 ~ know_covid*party_id + age_recode*party_id + gender_recode*party_id + race_recode*party_id + educ_recode*party_id + employment*party_id + income*party_id + insurance*party_id + inputstate*party_id + interview_date*party_id,
                      data     = subset(cces_model, year=="2020"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 2: Parties fully interacted
tr_pid_2 <- lm_robust(trump_2020 ~ covid_pc_int*party_id + age_recode*party_id + gender_recode*party_id + race_recode*party_id + educ_recode*party_id + employment*party_id + income*party_id + insurance*party_id + inputstate*party_id + interview_date*party_id,
                      data     = subset(cces_model, year=="2020"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 3: Parties fully interacted
tr_pid_3 <- lm_robust(trump_2020 ~ know_covid*party_id + covid_pc_int*party_id + age_recode*party_id + gender_recode*party_id + race_recode*party_id + educ_recode*party_id + employment*party_id + income*party_id + insurance*party_id + inputstate*party_id + interview_date*party_id,
                      data     = subset(cces_model, year=="2020"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1" = tr_pid_1,
                  "Model 2" = tr_pid_2,
                  "Model 3" = tr_pid_3),
             coef_omit   = "^(?!.*know_covid|.*covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = c("know_covid" =
                               "Respondent got COVID-19 or knows someone who did",
                             "know_covid:party_idRepublican" =
                               "Got COVID-19 or knows someone × Republican",
                             "party_idRepublican:know_covid" =
                               "Got COVID-19 or knows someone × Republican",
                             "know_covid:party_idDemocrat" =
                               "Got COVID-19 or knows someone × Democrat",
                             "party_idDemocrat:know_covid" =
                               "Got COVID-19 or knows someone × Democrat",
                             "know_covid:party_idOther/unsure" =
                               "Got COVID-19 or knows someone × Other/unsure",
                             "party_idOther/unsure:know_covid" =
                               "Got COVID-19 or knows someone × Other/unsure",
                             "covid_pc_int" =
                               "Cumulative county-level COVID-19 cases (5 per 100)",
                             "covid_pc_int:party_idRepublican" = 
                               "County-level cases (5 per 100) × Republican",
                             "party_idRepublican:covid_pc_int" = 
                               "County-level cases (5 per 100) × Republican",
                             "covid_pc_int:party_idDemocrat" = 
                               "County-level cases (5 per 100) × Democrat",
                             "party_idDemocrat:covid_pc_int" = 
                               "County-level cases (5 per 100) × Democrat",
                             "covid_pc_int:party_idOther/unsure" = 
                               "County-level cases (5 per 100) × Other/unsure",
                             "party_idOther/unsure:covid_pc_int" = 
                               "County-level cases (5 per 100) × Other/unsure"),
             add_rows    =
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Fully interacted Trump models for party identification.",
             stars       = FALSE,
             output      = "Trump table for parties interacted.tex") %>%
  row_spec(row=c(27,30), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Trump support" = 3))

##############################################################################
# Trump subgroup analyses: Strength of partisanship
##############################################################################

# Model 1: Strong partisans only
tr_str_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_20)),
                      data     = subset(cces_model, year=="2020" & party_lean=="Strong partisan"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 2: Strong partisans only
tr_str_2 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_20)),
                      data     = subset(cces_model, year=="2020" & party_lean=="Strong partisan"),
                      weights  = commonweight,
                      clusters = countyfips,
                      se_type  = "stata")

# Model 1: Weak/non-partisans only
tr_weak_1 <- lm_robust(as.formula(paste0("trump_2020 ~ know_covid", cov_20)),
                       data     = subset(cces_model, year=="2020" & party_lean=="Weak/non-partisan"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 2: Weak/non-partisans only
tr_weak_2 <- lm_robust(as.formula(paste0("trump_2020 ~ covid_pc_int", cov_20)),
                       data     = subset(cces_model, year=="2020" & party_lean=="Weak/non-partisan"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 1: Partisan strength fully interacted
tr_lean_1 <- lm_robust(trump_2020 ~ know_covid*party_lean + age_recode*party_lean + gender_recode*party_lean + race_recode*party_lean + educ_recode*party_lean + employment*party_lean + income*party_lean + insurance*party_lean + party_id*party_lean + inputstate*party_lean + interview_date*party_lean,
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 2: Partisan strength fully interacted
tr_lean_2 <- lm_robust(trump_2020 ~ covid_pc_int*party_lean + age_recode*party_lean + gender_recode*party_lean + race_recode*party_lean + educ_recode*party_lean + employment*party_lean + income*party_lean + insurance*party_lean + party_id*party_lean + inputstate*party_lean + interview_date*party_lean,
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Model 3: Partisan strength fully interacted
tr_lean_3 <- lm_robust(trump_2020 ~ know_covid*party_lean + covid_pc_int*party_lean + age_recode*party_lean + gender_recode*party_lean + race_recode*party_lean + educ_recode*party_lean + employment*party_lean + income*party_lean + insurance*party_lean + party_id*party_lean + inputstate*party_lean + interview_date*party_lean,
                       data     = subset(cces_model, year=="2020"),
                       weights  = commonweight,
                       clusters = countyfips,
                       se_type  = "stata")

# Compile results into table
modelsummary(list("Model 1" = tr_lean_1,
                  "Model 2" = tr_lean_2,
                  "Model 3" = tr_lean_3),
             coef_omit   = "^(?!.*know_covid|.*covid_pc_int)",
             fmt         = "%.3f",
             statistic   = c("({std.error})", "P={p.value}"),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_rename = c("know_covid" =
                               "Respondent got COVID-19 or knows someone who did",
                             "know_covid:party_leanStrong partisan" =
                               "Got COVID-19 or knows someone × Strong partisan",
                             "party_leanStrong partisan:know_covid" =
                               "Got COVID-19 or knows someone × Strong partisan",
                             "covid_pc_int" =
                               "Cumulative county-level COVID-19 cases (5 per 100)",
                             "covid_pc_int:party_leanStrong partisan" =
                               "County-level cases (5 per 100) × Strong partisan",
                             "party_leanStrong partisan:covid_pc_int" =
                               "County-level cases (5 per 100) × Strong partisan"),
             add_rows    =
               as.data.frame(
                 rbind(cbind("Controls",                     "Yes", "Yes", "Yes"),
                       cbind("State fixed effects",          "Yes", "Yes", "Yes"),
                       cbind("Interview date fixed effects", "Yes", "Yes", "Yes"),
                       cbind("Clustered standard errors", "County", "County", "County"),
                       cbind("Survey weights",               "Yes", "Yes", "Yes"))),
             title       = "Fully interacted Trump models for partisan strength.",
             stars       = FALSE,
             output      = "Trump table for strength interacted.tex") %>%
  row_spec(row=c(15,18), hline_after=TRUE) %>%
  add_header_above(c(" " = 1, "Trump support" = 3))

##############################################################################
# Graph: All subgroup analyses
##############################################################################

# Capture M4A subgroup coefficients
m4a_sub_vals <- data.frame(rbind(
  
  # Partisan identification
  cbind("know_covid", "Partisan identification", "Democrat",
        m4a_dem_1$coef["know_covid"], m4a_dem_1$std.error["know_covid"]),
  cbind("know_covid", "Partisan identification", "Republican",
        m4a_rep_1$coef["know_covid"], m4a_rep_1$std.error["know_covid"]),
  cbind("know_covid", "Partisan identification", "Independent",
        m4a_ind_1$coef["know_covid"], m4a_ind_1$std.error["know_covid"]),
  cbind("covid_pc_int", "Partisan identification", "Democrat",
        m4a_dem_2$coef["covid_pc_int"], m4a_dem_2$std.error["covid_pc_int"]),
  cbind("covid_pc_int", "Partisan identification", "Republican",
        m4a_rep_2$coef["covid_pc_int"], m4a_rep_2$std.error["covid_pc_int"]),
  cbind("covid_pc_int", "Partisan identification", "Independent",
        m4a_ind_2$coef["covid_pc_int"], m4a_ind_2$std.error["covid_pc_int"]),
  
  # Partisan strength
  cbind("know_covid", "Partisan strength", "Strong partisan",
        m4a_str_1$coef["know_covid"], m4a_str_1$std.error["know_covid"]),
  cbind("know_covid", "Partisan strength", "Weak/non-partisan",
        m4a_weak_1$coef["know_covid"], m4a_weak_1$std.error["know_covid"]),
  cbind("covid_pc_int", "Partisan strength", "Strong partisan",
        m4a_str_2$coef["covid_pc_int"], m4a_str_2$std.error["covid_pc_int"]),
  cbind("covid_pc_int", "Partisan strength", "Weak/non-partisan",
        m4a_weak_2$coef["covid_pc_int"], m4a_weak_2$std.error["covid_pc_int"])))
colnames(m4a_sub_vals) <- c("tx", "group", "level", "coef", "se")

# Reorder party identification
m4a_sub_vals$level <- factor(m4a_sub_vals$level, levels = c("Independent", "Republican", "Democrat", "Weak/non-partisan", "Strong partisan"))

# Treat values as numeric
m4a_sub_vals$coef <- as.numeric(m4a_sub_vals$coef)
m4a_sub_vals$se   <- as.numeric(m4a_sub_vals$se)

# Coefficient plots
m4a_sub_coef_plot_1 <- ggplot(subset(m4a_sub_vals, tx == "know_covid"),
                              aes(x=level, y=coef,
                                  ymin=coef-1.96*se, ymax=coef+1.96*se)) +
  geom_point() +
  geom_errorbar(width=0.2) +
  theme_test() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        panel.spacing = unit(1, "lines"),
        text = element_text(size = 10, face = "bold"),
        strip.background = element_blank(),
        strip.text = element_text(size = 10),
        axis.ticks   = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25)) +
  coord_flip() +
  ggtitle("") +
  ylab("\nPercentage-point difference in\nM4A support given getting COVID-19\nor knowing someone who did") +
  geom_hline(yintercept=0, color="black") +
  scale_y_continuous(limits=c(-0.05, 0.15), labels=scales::percent_format(accuracy=1)) +
  ggforce::facet_col(facets = vars(group), scales = "free_y", space = "free")

m4a_sub_coef_plot_2 <- ggplot(subset(m4a_sub_vals, tx == "covid_pc_int"),
                              aes(x=level, y=coef,
                                  ymin=coef-1.96*se, ymax=coef+1.96*se)) +
  geom_point() +
  geom_errorbar(width=0.2) +
  theme_test() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        panel.spacing = unit(1, "lines"),
        text = element_text(size = 10, face = "bold"),
        strip.background = element_blank(),
        strip.text = element_text(size = 10),
        axis.ticks   = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25)) +
  coord_flip() +
  ggtitle("") +
  ylab("\nPercentage-point difference in\nM4A support given 5 per 100 increase\nin COVID-19 cases per capita") +
  geom_hline(yintercept=0, color="black") +
  scale_y_continuous(limits=c(-0.05, 0.15), labels=scales::percent_format(accuracy=1)) +
  ggforce::facet_col(facets = vars(group), scales = "free_y", space = "free")

# Capture Trump subgroup coefficients
tr_sub_vals <- data.frame(rbind(
  
  # Partisan identification
  cbind("know_covid", "Partisan identification", "Democrat",
        tr_dem_1$coef["know_covid"], tr_dem_1$std.error["know_covid"]),
  cbind("know_covid", "Partisan identification", "Republican",
        tr_rep_1$coef["know_covid"], tr_rep_1$std.error["know_covid"]),
  cbind("know_covid", "Partisan identification", "Independent",
        tr_ind_1$coef["know_covid"], tr_ind_1$std.error["know_covid"]),
  cbind("covid_pc_int", "Partisan identification", "Democrat",
        tr_dem_2$coef["covid_pc_int"], tr_dem_2$std.error["covid_pc_int"]),
  cbind("covid_pc_int", "Partisan identification", "Republican",
        tr_rep_2$coef["covid_pc_int"], tr_rep_2$std.error["covid_pc_int"]),
  cbind("covid_pc_int", "Partisan identification", "Independent",
        tr_ind_2$coef["covid_pc_int"], tr_ind_2$std.error["covid_pc_int"]),
  
  # Partisan strength
  cbind("know_covid", "Partisan strength", "Strong partisan",
        tr_str_1$coef["know_covid"], tr_str_1$std.error["know_covid"]),
  cbind("know_covid", "Partisan strength", "Weak/non-partisan",
        tr_weak_1$coef["know_covid"], tr_weak_1$std.error["know_covid"]),
  cbind("covid_pc_int", "Partisan strength", "Strong partisan",
        tr_str_2$coef["covid_pc_int"], tr_str_2$std.error["covid_pc_int"]),
  cbind("covid_pc_int", "Partisan strength", "Weak/non-partisan",
        tr_weak_2$coef["covid_pc_int"], tr_weak_2$std.error["covid_pc_int"])))
colnames(tr_sub_vals) <- c("tx", "group", "level", "coef", "se")

# Reorder party identification
tr_sub_vals$level <- factor(tr_sub_vals$level, levels = c("Independent", "Republican", "Democrat", "Weak/non-partisan", "Strong partisan"))

# Treat values as numeric
tr_sub_vals$coef <- as.numeric(tr_sub_vals$coef)
tr_sub_vals$se   <- as.numeric(tr_sub_vals$se)

# Coefficient plots
tr_sub_coef_plot_1 <- ggplot(subset(tr_sub_vals, tx == "know_covid"),
                             aes(x=level, y=coef,
                                 ymin=coef-1.96*se, ymax=coef+1.96*se)) +
  geom_point() +
  geom_errorbar(width=0.2) +
  theme_test() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        panel.spacing = unit(1, "lines"),
        text = element_text(size = 10, face = "bold"),
        strip.background = element_blank(),
        strip.text = element_text(size = 10),
        axis.ticks   = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25)) +
  coord_flip() +
  ggtitle("") +
  ylab("\nPercentage-point difference in\nTrump support given getting COVID-19\nor knowing someone who did") +
  geom_hline(yintercept=0, color="black") +
  scale_y_continuous(limits=c(-0.15, 0.05), labels=scales::percent_format(accuracy=1)) +
  ggforce::facet_col(facets = vars(group), scales = "free_y", space = "free")

tr_sub_coef_plot_2 <- ggplot(subset(tr_sub_vals, tx == "covid_pc_int"),
                             aes(x=level, y=coef,
                                 ymin=coef-1.96*se, ymax=coef+1.96*se)) +
  geom_point() +
  geom_errorbar(width=0.2) +
  theme_test() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        panel.spacing = unit(1, "lines"),
        text = element_text(size = 10, face = "bold"),
        strip.background = element_blank(),
        strip.text = element_text(size = 10),
        axis.ticks   = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25)) +
  coord_flip() +
  ggtitle("") +
  ylab("\nPercentage-point difference in\nTrump support given 5 per 100 increase\nin COVID-19 cases per capita") +
  geom_hline(yintercept=0, color="black") +
  scale_y_continuous(limits=c(-0.15, 0.05), labels=scales::percent_format(accuracy=1)) +
  ggforce::facet_col(facets = vars(group), scales = "free_y", space = "free")

# Combine plots together
sub_plots_all <-
  cowplot::plot_grid(m4a_sub_coef_plot_1, m4a_sub_coef_plot_2,
                     tr_sub_coef_plot_1,  tr_sub_coef_plot_2,
                     labels=c("A. Medicare for All support", "",
                              "B. Trump support", ""), hjust=0, label_x=0.05)

# Print figure
ggsave(plot=sub_plots_all, file="Subgroup coefficient plots.png",
       width=7.5, height=6.5, units='in', dpi=600)

##############################################################################
# Graph: Medicare for All support by state
##############################################################################

# Support by state
states_2019 <- subset(cces_model, year=="2019") %>%
  group_by(state_alt) %>%
  summarise(mean = weighted.mean(medi_exp, commonweight, na.rm=T))
states_2019$state <- states_2019$state_alt

states_2019 <- states_2019 %>% mutate(
  cats = case_when(
    mean < 0.40 ~ "Less than 40%",
    mean < 0.50 ~ "40% to 49%",
    mean < 0.60 ~ "50% to 59%",
    mean < 0.70 ~ "60% to 69%",
    mean < 1    ~ "70% or higher",
  ))
states_2019$cats <- factor(states_2019$cats, levels = c("70% or higher", "60% to 69%", "50% to 59%", "40% to 49%", "Less than 40%"))

states_2020 <- subset(cces_model, year=="2020") %>%
  group_by(inputstate) %>%
  summarise(mean = weighted.mean(medi_exp, commonweight, na.rm=T))
states_2020$fips <- states_2020$inputstate

states_2020 <- states_2020 %>% mutate(
  cats = case_when(
    mean < 0.40 ~ "Less than 40%",
    mean < 0.50 ~ "40% to 49%",
    mean < 0.60 ~ "50% to 59%",
    mean < 0.70 ~ "60% to 69%",
    mean < 1    ~ "70% or higher",
  ))
states_2020$cats <- factor(states_2020$cats, levels = c("70% or higher", "60% to 69%", "50% to 59%", "40% to 49%", "Less than 40%"))

# Map for 2019
map_medi_19 <- plot_usmap(regions="states", data=states_2019,
                          values="cats", size=0.4) +
  ggtitle("Late 2019") +
  theme(legend.position = "none",
        plot.title = element_text(size=14, face="bold", hjust=0.5)) +
  scale_fill_manual(values=c("#440154FF", "#39568CFF", "#1F968BFF", "#73D055FF", "#FDE725FF"))

# Map for 2020
map_medi_20 <- plot_usmap(regions="states", data=states_2020,
                          values="cats", size=0.4) +
  ggtitle("Late 2020") +
  theme(legend.position = "none",
        plot.title = element_text(size=14, face="bold", hjust=0.5)) +
  scale_fill_manual(values=c("#440154FF", "#39568CFF", "#1F968BFF", "#73D055FF", "#FDE725FF"))

# Legend
legend <- plot_usmap(regions="states", data=states_2019, values="cats") +
  theme(legend.position = "right",
        legend.text = element_text(size = 9),
        legend.title = element_text(size = 10, face="bold"),
        legend.background = element_blank()) +
  scale_fill_manual(name="Support for\nMedicare for All",
                    values=c("#440154FF", "#39568CFF", "#1F968BFF", "#73D055FF", "#FDE725FF"))

# Compile figures into object
medicare_maps <- plot_grid(
  map_medi_19, map_medi_20,
  plot_grid(get_legend(legend), NULL, nrow=2, rel_heights = c(1,0.2)),
  nrow=1, rel_widths=c(1,1,0.5))

# Print figure
ggsave(plot=medicare_maps, file="Medicare for All map.png",
       width=7, height=2.5, units='in', dpi=600)

##############################################################################
# Graph: COVID-19 case rates by county
##############################################################################

# Calculate cases per capita
# Choose date near median interview date (±6 months)
covid$DP05_0001E <- as.numeric(as.character(covid$DP05_0001E))
covid$case_rate_apr_15_20 <- covid$`2020-04-15`/covid$DP05_0001E
covid$case_rate_oct_15_20 <- covid$`2020-10-15`/covid$DP05_0001E
covid$case_rate_apr_15_21 <- covid$`2021-04-15`/covid$DP05_0001E

# Select relevant columns
covid_subset <- covid %>% select(fips, case_rate_apr_15_20, case_rate_oct_15_20, case_rate_apr_15_21)

# Map for April 15, 2020
map_apr_15_20 <- plot_usmap(regions="counties", data=covid_subset,
                            values="case_rate_apr_15_20", color="NA") +
  ggtitle("April 15, 2020") +
  theme(legend.position = "none",
        plot.title = element_text(size=14, face="bold", hjust=0.5)) +
  scale_fill_viridis(direction=-1, limits=c(0,0.2))

# Map for October 15, 2020
map_oct_15_20 <- plot_usmap(regions="counties", data=covid_subset,
                            values="case_rate_oct_15_20", color="NA") +
  ggtitle("October 15, 2020") +
  theme(legend.position = "none",
        plot.title = element_text(size=14, face="bold", hjust=0.5)) +
  scale_fill_viridis(direction=-1, limits=c(0,0.2), oob=squish)

# Map for April 15, 2021
map_apr_15_21 <- plot_usmap(regions="counties", data=covid_subset,
                            values="case_rate_apr_15_21", color="NA") +
  ggtitle("April 15, 2021") +
  theme(legend.position = "none",
        plot.title = element_text(size=14, face="bold", hjust=0.5)) +
  scale_fill_viridis(direction=-1, limits=c(0,0.2), oob=squish)

# Legend
legend <- plot_usmap(regions="counties", data=covid_subset,
                     values="case_rate_apr_15_21") +
  theme(legend.position = "right",
        legend.text = element_text(size = 9),
        legend.title = element_text(size = 10, face="bold"),
        legend.background = element_blank()) +
  scale_fill_viridis(direction=-1, limits=c(0,0.15), breaks=c(0,0.05,0.1,0.15), labels=c("  0 cases"," 5 per 100", "10 per 100", "15 per 100\n or above"), oob=squish, name="COVID-19 case rate\n")

# Compile figures into object
covid_maps <- plot_grid(
  map_apr_15_20,
  plot_grid(get_legend(legend), NULL, nrow=2, rel_heights=c(1,0.22)),
  map_oct_15_20,
  plot_grid(get_legend(legend), NULL, nrow=2, rel_heights=c(1,0.22)),
  map_apr_15_21,
  plot_grid(get_legend(legend), NULL, nrow=2, rel_heights=c(1,0.22)),
  nrow=3, rel_widths=c(1,0.5))

# Print figure
ggsave(plot=covid_maps, file="COVID-19 cases map.png",
       width=5, height=8, units='in', dpi=600)
